This script shows how to analyze a typical ISS segmented dataset. The analysis is based on Scanpy (documentation in https://scanpy.readthedocs.io/en/stable/index.html) and Squidpy (documentation in https://squidpy.readthedocs.io/en/stable/)
import ISS_postprocessing
import scanpy as sc
import numpy as np
import matplotlib.pyplot as plt
>>>> using CPU Running test snippet to check if MKL-DNN working see https://pytorch.org/docs/stable/backends.html?highlight=mkl ** MKL version working - CPU version is sped up. ** >>>> using CPU Running test snippet to check if MKL-DNN working see https://pytorch.org/docs/stable/backends.html?highlight=mkl ** MKL version working - CPU version is sped up. ** DIPlib -- a quantitative image analysis library Version 3.2.0 (Feb 8 2022) For more information see https://diplib.org
samples = [
'/media/christoffer/hfsc_processing_2/HFSC_R5/',
'/media/christoffer/hfsc_processing_2/HFSC_R6/',
'/media/christoffer/hfsc_processing_2/HFSC_R7/',
'/media/christoffer/hfsc_processing_2/HFSC_R8/'
]
for sample in samples:
ISS_postprocessing.segmentation.stardist_segmentation(image_path = sample + '/preprocessing/stitched/Round4_4.tif',
output_path = sample,
model_name = '2D_versatile_fluo',
expand_cells = True,
n_tiles = (4,4),
expanded_distance = 20,)
2022-06-14 16:15:37.294989: I tensorflow/core/util/util.cc:169] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`. 2022-06-14 16:15:37.298974: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /home/christoffer/anaconda3/envs/3_ISS_postprocessing_new/lib/python3.8/site-packages/cv2/../../lib64: 2022-06-14 16:15:37.298992: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.
Found model '2D_versatile_fluo' for 'StarDist2D'. Loading network weights from 'weights_best.h5'.
2022-06-14 16:15:38.588043: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcuda.so.1'; dlerror: libcuda.so.1: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /home/christoffer/anaconda3/envs/3_ISS_postprocessing_new/lib/python3.8/site-packages/cv2/../../lib64: 2022-06-14 16:15:38.588067: W tensorflow/stream_executor/cuda/cuda_driver.cc:269] failed call to cuInit: UNKNOWN ERROR (303) 2022-06-14 16:15:38.588097: I tensorflow/stream_executor/cuda/cuda_diagnostics.cc:156] kernel driver does not appear to be running on this host (borka): /proc/driver/nvidia/version does not exist 2022-06-14 16:15:38.588323: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations: AVX2 AVX512F AVX512_VNNI FMA To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.
Loading thresholds from 'thresholds.json'. Using default values: prob_thresh=0.479071, nms_thresh=0.3. (7701, 28291) normalize image predict instance 1/1 [==============================] - 0s 319ms/step 1/1 [==============================] - 0s 82ms/step
0%| | 0/16 [00:00<?, ?it/s]
1/1 [==============================] - 6s 6s/step
6%|███████████▏ | 1/16 [00:06<01:32, 6.17s/it]
1/1 [==============================] - 4s 4s/step
12%|██████████████████████▎ | 2/16 [00:11<01:16, 5.45s/it]
1/1 [==============================] - 5s 5s/step
19%|█████████████████████████████████▍ | 3/16 [00:16<01:10, 5.46s/it]
1/1 [==============================] - 5s 5s/step
25%|████████████████████████████████████████████▌ | 4/16 [00:21<01:04, 5.39s/it]
1/1 [==============================] - 5s 5s/step
31%|███████████████████████████████████████████████████████▋ | 5/16 [00:27<00:59, 5.39s/it]
1/1 [==============================] - 5s 5s/step
38%|██████████████████████████████████████████████████████████████████▊ | 6/16 [00:32<00:54, 5.41s/it]
1/1 [==============================] - 5s 5s/step
44%|█████████████████████████████████████████████████████████████████████████████▉ | 7/16 [00:37<00:47, 5.30s/it]
1/1 [==============================] - 5s 5s/step
50%|█████████████████████████████████████████████████████████████████████████████████████████ | 8/16 [00:43<00:42, 5.32s/it]
1/1 [==============================] - 5s 5s/step
56%|████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 9/16 [00:48<00:37, 5.31s/it]
1/1 [==============================] - 5s 5s/step
62%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 10/16 [00:54<00:33, 5.51s/it]
1/1 [==============================] - 5s 5s/step
69%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 11/16 [00:59<00:27, 5.47s/it]
1/1 [==============================] - 5s 5s/step
75%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 12/16 [01:04<00:21, 5.35s/it]
1/1 [==============================] - 5s 5s/step
81%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 13/16 [01:10<00:16, 5.39s/it]
1/1 [==============================] - 5s 5s/step
88%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 14/16 [01:16<00:11, 5.56s/it]
1/1 [==============================] - 5s 5s/step
94%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 15/16 [01:21<00:05, 5.55s/it]
1/1 [==============================] - 4s 4s/step
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 16/16 [01:26<00:00, 5.41s/it]
label image expand image save output Found model '2D_versatile_fluo' for 'StarDist2D'. Loading network weights from 'weights_best.h5'. Loading thresholds from 'thresholds.json'. Using default values: prob_thresh=0.479071, nms_thresh=0.3. (7700, 28269) normalize image predict instance 1/1 [==============================] - 0s 191ms/step 1/1 [==============================] - 0s 78ms/step
0%| | 0/16 [00:00<?, ?it/s]
1/1 [==============================] - 6s 6s/step
6%|███████████▏ | 1/16 [00:06<01:30, 6.05s/it]
1/1 [==============================] - 5s 5s/step
12%|██████████████████████▎ | 2/16 [00:11<01:18, 5.64s/it]
1/1 [==============================] - 4s 4s/step
19%|█████████████████████████████████▍ | 3/16 [00:16<01:07, 5.20s/it]
1/1 [==============================] - 5s 5s/step
25%|████████████████████████████████████████████▌ | 4/16 [00:21<01:03, 5.31s/it]
1/1 [==============================] - 5s 5s/step
31%|███████████████████████████████████████████████████████▋ | 5/16 [00:27<01:00, 5.51s/it]
1/1 [==============================] - 5s 5s/step
38%|██████████████████████████████████████████████████████████████████▊ | 6/16 [00:33<00:55, 5.54s/it]
1/1 [==============================] - 5s 5s/step
44%|█████████████████████████████████████████████████████████████████████████████▉ | 7/16 [00:38<00:50, 5.60s/it]
1/1 [==============================] - 5s 5s/step
50%|█████████████████████████████████████████████████████████████████████████████████████████ | 8/16 [00:44<00:45, 5.66s/it]
1/1 [==============================] - 6s 6s/step
56%|████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 9/16 [00:50<00:40, 5.85s/it]
1/1 [==============================] - 5s 5s/step
62%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 10/16 [00:56<00:35, 5.89s/it]
1/1 [==============================] - 6s 6s/step
69%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 11/16 [01:02<00:29, 5.96s/it]
1/1 [==============================] - 5s 5s/step
75%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 12/16 [01:08<00:23, 5.80s/it]
1/1 [==============================] - 6s 6s/step
81%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 13/16 [01:14<00:17, 5.87s/it]
1/1 [==============================] - 6s 6s/step
88%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 14/16 [01:20<00:11, 5.94s/it]
1/1 [==============================] - 5s 5s/step
94%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 15/16 [01:26<00:05, 5.94s/it]
1/1 [==============================] - 5s 5s/step
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 16/16 [01:31<00:00, 5.74s/it]
label image expand image save output Found model '2D_versatile_fluo' for 'StarDist2D'. Loading network weights from 'weights_best.h5'. Loading thresholds from 'thresholds.json'. Using default values: prob_thresh=0.479071, nms_thresh=0.3. (9564, 20780) normalize image predict instance 1/1 [==============================] - 0s 198ms/step 1/1 [==============================] - 0s 70ms/step
0%| | 0/16 [00:00<?, ?it/s]
1/1 [==============================] - 5s 5s/step
6%|███████████▏ | 1/16 [00:05<01:18, 5.24s/it]
1/1 [==============================] - 5s 5s/step
12%|██████████████████████▎ | 2/16 [00:10<01:15, 5.38s/it]
1/1 [==============================] - 5s 5s/step
19%|█████████████████████████████████▍ | 3/16 [00:15<01:09, 5.32s/it]
1/1 [==============================] - 5s 5s/step
25%|████████████████████████████████████████████▌ | 4/16 [00:21<01:02, 5.24s/it]
1/1 [==============================] - 4s 4s/step
31%|███████████████████████████████████████████████████████▋ | 5/16 [00:25<00:56, 5.10s/it]
1/1 [==============================] - 5s 5s/step
38%|██████████████████████████████████████████████████████████████████▊ | 6/16 [00:31<00:52, 5.25s/it]
1/1 [==============================] - 5s 5s/step
44%|█████████████████████████████████████████████████████████████████████████████▉ | 7/16 [00:36<00:47, 5.31s/it]
1/1 [==============================] - 5s 5s/step
50%|█████████████████████████████████████████████████████████████████████████████████████████ | 8/16 [00:42<00:43, 5.47s/it]
1/1 [==============================] - 5s 5s/step
56%|████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 9/16 [00:48<00:38, 5.45s/it]
1/1 [==============================] - 4s 4s/step
62%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 10/16 [00:52<00:31, 5.26s/it]
1/1 [==============================] - 5s 5s/step
69%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 11/16 [00:58<00:27, 5.45s/it]
1/1 [==============================] - 5s 5s/step
75%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 12/16 [01:04<00:21, 5.47s/it]
1/1 [==============================] - 5s 5s/step
81%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 13/16 [01:09<00:16, 5.40s/it]
1/1 [==============================] - 5s 5s/step
88%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 14/16 [01:15<00:10, 5.46s/it]
1/1 [==============================] - 5s 5s/step
94%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 15/16 [01:20<00:05, 5.53s/it]
1/1 [==============================] - 5s 5s/step
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 16/16 [01:26<00:00, 5.40s/it]
label image expand image save output Found model '2D_versatile_fluo' for 'StarDist2D'. Loading network weights from 'weights_best.h5'. Loading thresholds from 'thresholds.json'. Using default values: prob_thresh=0.479071, nms_thresh=0.3. (9565, 20780) normalize image predict instance 1/1 [==============================] - 0s 192ms/step 1/1 [==============================] - 0s 78ms/step
0%| | 0/16 [00:00<?, ?it/s]
1/1 [==============================] - 6s 6s/step
6%|███████████▏ | 1/16 [00:06<01:32, 6.18s/it]
1/1 [==============================] - 5s 5s/step
12%|██████████████████████▎ | 2/16 [00:12<01:24, 6.01s/it]
1/1 [==============================] - 4s 4s/step
19%|█████████████████████████████████▍ | 3/16 [00:16<01:09, 5.38s/it]
1/1 [==============================] - 5s 5s/step
25%|████████████████████████████████████████████▌ | 4/16 [00:21<01:03, 5.33s/it]
1/1 [==============================] - 6s 6s/step
31%|███████████████████████████████████████████████████████▋ | 5/16 [00:28<01:01, 5.63s/it]
1/1 [==============================] - 5s 5s/step
38%|██████████████████████████████████████████████████████████████████▊ | 6/16 [00:33<00:56, 5.65s/it]
1/1 [==============================] - 5s 5s/step
44%|█████████████████████████████████████████████████████████████████████████████▉ | 7/16 [00:39<00:51, 5.67s/it]
1/1 [==============================] - 5s 5s/step
50%|█████████████████████████████████████████████████████████████████████████████████████████ | 8/16 [00:44<00:44, 5.56s/it]
1/1 [==============================] - 4s 4s/step
56%|████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 9/16 [00:49<00:37, 5.31s/it]
1/1 [==============================] - 5s 5s/step
62%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 10/16 [00:54<00:31, 5.29s/it]
1/1 [==============================] - 5s 5s/step
69%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 11/16 [01:00<00:27, 5.44s/it]
1/1 [==============================] - 5s 5s/step
75%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 12/16 [01:06<00:21, 5.49s/it]
1/1 [==============================] - 5s 5s/step
81%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 13/16 [01:11<00:16, 5.49s/it]
1/1 [==============================] - 5s 5s/step
88%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 14/16 [01:17<00:11, 5.54s/it]
1/1 [==============================] - 5s 5s/step
94%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 15/16 [01:22<00:05, 5.50s/it]
1/1 [==============================] - 4s 4s/step
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 16/16 [01:27<00:00, 5.48s/it]
label image expand image save output
for sample_path in samples:
ad = ISS_postprocessing.annotated_objects.create_anndata_obj(spots_file = sample_path+'/decoded.csv',
segmentation_mask = sample_path+'stardist_segmentation_expanded.npz',#'cell_segmentation/cellpose_segmentation_expanded_2.npz'
output_file = sample_path+'/anndata_stardist.h5ad',
filter_data=False,
metric = 'distance',
write_h5ad = True,
value= 0.4,
convert_coords = True,
conversion_factor = 1)
reading spots file load coo file assign spots to cells
AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([2, 3], dtype='int64'), … AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([2, 3], dtype='int64'), …
write h5ad reading spots file load coo file
AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([0, 1], dtype='int64'), … AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([0, 1], dtype='int64'), …
assign spots to cells write h5ad reading spots file load coo file
AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([1, 2], dtype='int64'), … AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([1, 2], dtype='int64'), …
assign spots to cells write h5ad reading spots file load coo file
AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([0, 1], dtype='int64'), … AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([0, 1], dtype='int64'), …
assign spots to cells write h5ad
ad = ISS_postprocessing.annotated_objects.concat_anndata(sample_anndata_list = samples,
anndata_name = 'anndata_stardist.h5ad'
)
ad.obs['sample_id'] = ad.obs['sample_id'].str.split('/', expand = True)[4]
/media/christoffer/hfsc_processing_2/HFSC_R5/anndata_stardist.h5ad /media/christoffer/hfsc_processing_2/HFSC_R6/anndata_stardist.h5ad /media/christoffer/hfsc_processing_2/HFSC_R7/anndata_stardist.h5ad /media/christoffer/hfsc_processing_2/HFSC_R8/anndata_stardist.h5ad
Observation names are not unique. To make them unique, call `.obs_names_make_unique`.
sc.pp.calculate_qc_metrics(ad, percent_top=None, log1p=False, inplace=True)
keep = (ad.obs['total_counts'] > 5) & (ad.obs['n_genes_by_counts'] > 3)
print(sum(keep))
75168
ad = ad[keep,:]
We now save the raw data in adata.raw and apply normalization and log-transformation to minimize the effect of outliers. We also scale the data to give the same importance to all genes in clustering
</font>
ad.layers["raw"] = ad.X.copy()
sc.experimental.pp.normalize_pearson_residuals(ad)
normalized_exp = np.array(ad.X)
normalized_exp[normalized_exp< 0] = 0
ad.X = normalized_exp
sc.pp.log1p(ad)
Observation names are not unique. To make them unique, call `.obs_names_make_unique`. Observation names are not unique. To make them unique, call `.obs_names_make_unique`.
sc.tl.pca(ad, svd_solver='arpack')
plt.rcParams['figure.facecolor'] = 'white'
sc.pl.pca_variance_ratio(ad, log=True)
sc.pp.neighbors(ad, n_neighbors=20, n_pcs=40)
sc.tl.umap(ad,min_dist=0.005)
for i in [2, 3]:
if "norm_leiden_"+str(i) in ad.obs.columns:
plt.rcdefaults()
with plt.rc_context({'figure.figsize': (15, 10)}):
sc.pl.umap(ad, color = ("norm_leiden_"+str(i)),s=5,add_outline=True,legend_loc='on data',legend_fontsize=20,legend_fontoutline=2,alpha = 1)
else:
print('clustering @ resolution ' + str(i))
sc.tl.leiden(ad, resolution =i, key_added = ("norm_leiden_"+str(i)))
plt.rcdefaults()
with plt.rc_context({'figure.figsize': (15, 10)}):
sc.pl.umap(ad, color = ("norm_leiden_"+str(i)),s=5,add_outline=True,legend_loc='on data',legend_fontsize=20,legend_fontoutline=2)
spatial = np.array(ad.obs[['x','y']].astype('<f8'))
ad.obsm['spatial'] = spatial
ad.obsm['xy_loc'] = spatial
cluster = 'norm_leiden_2'
for sample in ad.obs['sample_id'].unique():
print(sample)
ad_int = ad[ad.obs['sample_id'] == sample]
with plt.rc_context({'figure.figsize': (25, 17.5)}):
sc.pl.spatial(ad_int,color=cluster ,spot_size=50)
HFSC_R5
HFSC_R6
HFSC_R7
HFSC_R8
sc.tl.rank_genes_groups(ad, cluster, method='wilcoxon', key_added = "wilcoxon")
sc.pl.rank_genes_groups(ad, n_genes=30, sharey=False, key="wilcoxon")
plt.rcdefaults()
sc.pl.rank_genes_groups_matrixplot(ad, n_genes=3, key='wilcoxon', groupby=cluster, cmap = 'viridis', )
ad.obs['project'] = 'sc'
for cluster_ in sorted(ad.obs[cluster].value_counts().head(5).index.astype(int)):
print(cluster_)
genes = list(sc.get.rank_genes_groups_df(ad,group=str(cluster_), key='wilcoxon')['names'].head(4))
genes.insert(0, cluster)
genes.insert(1, 'sample_id')
plt.rcdefaults()
sc.pl.umap(ad, color = genes, cmap = 'turbo', ncols = 3, legend_loc='on data',legend_fontsize=10,legend_fontoutline=2)
ISS_postprocessing.annotated_objects.plot_specific_cluster(ad,
clusters_to_map = cluster,
cluster = str(cluster_),
broad_cluster = 'project',
key='wilcoxon',
size=0.5,
number_of_marker_genes=10,
sample_id_column='sample_id',
dim_subplots=[4, 4],)
plt.show()
0
PAX3 EGFLAM RSPO3 CCL21 SPRY1 SIM1 FGF9 ID4 NKX2-8 SPON1
1
HOXA9 GBX2 SIM1 PTPRZ1 VEGFC NOTCH2 FGF20 COL12A1 OLIG2 HOXA11
2
NKX6-2 SPON1 NOTCH2 COL12A1 PTPRZ1 WNT3A FGF9 PDGFD FGF20 NKX2-8
3
FGF7 SPON1 WNT3A ID4 EGFLAM HOXA2 RSPO3 NKX2-8 FGF9 NOTCH2
4
GSX2 RSPO3 EGFLAM NOTCH2 WNT3A ID4 SPON1 SIM1 HOXA2 FGF9